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ABSTRACT 


The complex LMS adaptive algorithm developed by Widrow 
[Ref. 1] is used in the frequency domain to estimate the 
azimuth and elevation angles of a plane wave incident upon 
a planar array. The complex LMS algorithm is applied to two 
cases. The first case is a passive detection problem. The 
second case is a pulse communication problem. [In both 
cases, complex weights are determined using the complex LMS 
Pesonlehnm whtenecopnase all of the output electrical signals 
from.the planar array. Three versions of the complex LMS 


algorithm are studied and their performances are compared. 
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Frequency domain beamforming is accomplished by applying 
appropriate phase shifts at the sensor outputs of an array 
Pemaccount for the selative propagation delays ofa signal 
from a particular direction. The phase-shifted signals from 
all sensors are then added together coherently to realize 
the full array gain. Discrete Fourier Transform (DFT) beam- 
forming is the usual method of determining the direction of 
arrival of a plane wave signal. A discrete number of 
spatial frequency bins are formed and each bin corresponds 
to a discrete direction. If the number of spatial frequency 
bins is large, very fine spatial resolution can be obtained. 

The phase shifts needed to cancel the relative pEOpaAgas 
tion delays can be determined adaptively. The complex LMS 
adaptive algorithm is used in this thesis. The LMS adaptive 
filter adjusts its adaptive weights recursively to minimize 
the mean square difference between a reference signal and 
its estimate. When the beam is steered toward a signal 
pEeMpogamimnage in a particular GQirection, the phase wf the 
Signals at all sensors must be the same. Therefore, the 
Signal at any sensor can be used as a reference which the 
others will be matched. The estimated signal is obtained by 
weighting the input signal by the current adaptive weights. 


Note that no prior knowledge of the reference signal is 


1: 


required. The response of the LMS adaptive filter converges 
to the discrete Wiener filter without a priori knowledge of 
the input [Refs. 1,2]. [Ref. 1] proposed the complex LMS 
algorithm to deal with complex inputs. (Ref. 5] addressed 
the implementation of the LMS adaptive filter in the fre- 
quency domain. [Ref. 4] used the LMS adaptive filter in the 
frequency domain to estamate the bearang of sauplanemwae 

due to an acoustic source radiating a sinusoidal signal. Buia 
this application, the LMS adaptive filter was implemented to 
estimate the phase difference between two sonar arrays 
separated by a distance many times the signal wavelength. 

The angle of arrival of a plane wave can be estimated if the 
frequency of the acoustic signal and the speed of wave propa- 
gation are known or can be extracted from the received signal 
itself. 

The objective of this thesis is to extend the results in 
[Ref. 3] and [Ref. 4] to a planar array of M xN elements 
(hydrophones) where M and N are greater than two. Such an 
array has an overall size many times the wavelength of the 
received signal. The inter-element separation, however, is 
usually maintained at a distance of less than or equal to 
one-half of the expected minimum wavelength. This requirement 
[Ref. 5] prevents the occurrence of grating lobes in the 
far-field beam pattern. A two-dimensional array allows 
spatial resolution in beth azimuth and elevation lei. a 


though the detection range in underwater acoustics is large 


eZ 


compared to the ocean depth, the effect of ray bending due 
to the inhomogeneous ocean medium can bend the incident 
acoustic wave such that it can appear to arrive at a steeper 
or shallower angle than the line of sight angle in the homo- 
geneous medium case. The elevation/depression angle is 
at present used to estimate Convergence Zone (CZ) and Bottom 
Bowumce (BB) ranges. 

in this thesis, the problem of a — wave incident upon 
a planar array of M xN elements is studied. The acoustic 
wave signal at each of the elements in the array are identical 
if the array is steered in the direction of the incident 
wavefront. However, if the main lobe of the array is not 
steered properly, the plane wave signal will have the same 
Spectral content at each element but modified by a phase 
Dibhee pisepemetonalmre: the slocation of the element with respect 
to some reference element. These undesirable phase shifts 
can be cancelled by applying appropriate phase weights at 
each element and thereby cophasing the total array output to 
realize its array gain. [Ref. 4] demonstrated that the LMS 
adaptive filter can achieve phase alignment between a refer- 
ence signal and input signal in the frequency domain by 
deeeeedpolicatcton of the complex LMS algorithm. This is 
equivalent to a tapped delay line structure in the time do- 
main. However, in the frequency domain, a time delay T 
corresponds to multiplication of a complex number that is 
ga 


equal to e where w is the signal frequency. The 


AS 


implementation of the LMS adaptive filter in the frequency 
domain requires fewer computations per iteration than in the 
time domain. An added advantage of using phase weighting 

is that a continuous range of spatial directions can be 
described, whereas in a tapped delay line structure, only 
finite increments of delays can be applied. Figure 1 shows 
a functional diagram of an N-element adaptive array imple- 
mented in the frequency domain [Ref. 6]. 

Chapter II of this thesis describes the specific struc- 
ture of the adaptive filter and the equations implemented for 
Simulations. The assumptions made in.the model are discussed 
and justified. Several modifications to the complex LMS 
adaptive algorithms are made to increase the array's spatial 
coverage and to ensure that the steady state phase weights 
do correspond to the direction of the incident wave. 

In Chapter III, a passive sonar system is modeled to 
test the ability of the modified complex LMS algorithm to 
estimate the direction of a source in the presence 6f noises 
The simulation program is implemented in VS APL. 

Chapter IV demonstrates the application of the complex 
LMS algorithm to a pulse communication problem. The inte- 
gration time in this case 1s much shorter than that of 3am 
passive sonar case. Two types of pulse waveforms are included, 
continuous wave (CW) and linear frequency-modulated (LFM). 


This simulation program is implemented in VS FORTRAN. 
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Figure I. Frequency Domain Adaptive Array 
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Chapter V concludes this thesis by identifying further 
research areas and other possible applications. 

Appendix A contains the derivation of the complex LMS 
algorithm. Appendix B has the description of the passive 
detection program implemented in APL. Appendix C has the 
description of the pulse communication program implemented 


in VS Fortran. 
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Po OVERVIEW OF ARRAYS (Refs. 5,6,7,8] 

The characteristics of the array elements and their 
arrangement in forming the array determine the ultimate 
P-GLormance GL an adaptive array system [Ref. 6]. Both the 
linear array and the planar array are examined here. 

1. Linear Arrays tRefs. Sy on 

Consider a linear array that has M equally spaced, 
identical point source elements along the x-axis. For illus- 
tration, Figure 2 shows a 7-element linear array with uniform 
interelement spacing d,, and a plane wave arriving at the 


array with an incident angle 9 as measured from the array 


normal. The phasor sum of all elements is: 

M~-1 a 

sm = } yit)e!™ (2.1) 
m=0 
, , ea th 

where wW is the phase shift between the m and the m+l 
elemett Seetem = O7,d2535, ...p0M-1. 

d.. sing aasind 

_ X ey 2s 

where: 


alee, 





Figure 2. 


mak mss meg mes mre 


A Seven-Element Linear Array with Uniform 
Interelement Spacing and Point Source Elements 
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k = wave number in radians per meter 


X = wavelength in meters 
f = frequency in hertz 
d,. = interelement spacing in meters 


= incident angle measured from array normal 


c = speed of wave propagation in meters per second. 
Rewriting equation (2.2) by substituting = c/f yields 


Gives Sine 





C C 


ie HOurier transform Of equation (2.1) is: 


S(f) = Fis(t)} = f ees.) at (2.4) 
co M-] ; 
= ff J yltyedMe7IeTE tay C5) 
-co m=0 
M~1 00 
= jy @Jmy f loa JOU ERs (2.6) 
mo 00 
S(f£,9) = Alt,8)Y(f£) (rou) 


where Y(f) is the frequency spectrum of the incident wave and 


A(f,8) is called the space factor or array factor. The array 


ne, 


factor A(f,@) determines the directional plaw] Clee iesaeea, 
in a plane containing the array. The dependence of A(f,@) on 
frequency, speed of propagation, element separation, number 


of elements, and incident angle can be shown by rewriting 


A(£,8) as: 
“eZ 
M-l, m-1 jm “5 d,sin 6) 
A(f,6) = jy) et™ = Je °S (aes 
m= 0 i 
Summing equation (2.8) yields: 
M 
Iv Sin >V 
A(f,8) = e ey (2295) 
Gana 7 


The normalized directlonal pattern vis semen. 


2 | 
ele 9) = ie log , {ates 8) } (2 a 


Me 
For nonisotropic (non-point source) elements, it 1s necessary 
to introduce an additional factor E(f,9) in equation (2aee 
to include the directional response pattern introduced by 
each sensor element [Ref. 6]. The overall directivity pat- 
tern then is given by the produce of the array factor and the 


element factor [Ref. 5]. However, if the size of the individual 
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elements are small compared to a wavelength, they can be 
assumed to be omnidirectional point sources, i.e., E(f,98) =1. 
The effects of increasing the number of elements while main- 
taining the same element spacing are shown in Figures 3 and 
4 [Ref. 6]. It can be seen that the main lob beamwidth de- 
Greases as the number of elements increases, and the number 
Of sidelobes and nulls increases. In Figures 5-8, the number 
of elements is held constant at M = 7 while the interelement 
Spacing 1S varied to illustrate the effects of elements 
Spacing on the directivity pattern. 

Beamsteering [Ref. 5] is accomplished by applying a 
linear phase shift across the line array as shown in Figure 
9 [Ref. 6]. The effect of the insertion of this sequence of 
phase shifts is that the main lobe is steered to an angle as 


measured off the boresight equal to 6 where 


g = sin *{i 5} (2.11) 
2d 
Xx 
and 6 is the phase shift between adjacent elements. Figure 


10 shows the directivity pattern of a steered linear array. 
7 onameaerdys TREES. 5,6) 
Much of the analysis done in linear arrays can be 
extended to the case of a rectangular-shaped planar array. 
A circular planar array or a spherical volume array would re- 
quire the use of polar and spherical coordinates respectively. 
However, array theory is invariant under coordinate 


transformations. 
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Figure 3. Directivity Pattern of a Three-Element 
Linear Array [Ref. 6:p. 40] 
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Figure 4. Directivity eateermnsetec Our EL emenia 
Linear Array [Ref. 6:p. 40] 
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Figure 5. Directivity Pattern of a Seven-Element 
Linear Array for d/A = 0.1 [Ref. 6:p. 41] 
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PiMPeem ome necctlylty Pattern Of a Seven-Element 
Linear Array for d/A = 0.2 [{Ref. 6:p. 41] 
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Figure 7. Directivity Pattern of a Seven-Element 


Linear Array for d/i = 0.5 [{Ref. 6:p. 42] 
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Figure 8. Directivity Pattern of a Seven-Element 
Linear Array for d =  [Ref. 6:p. 42] 
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Figure 9. Beamsteering by Applying a Linear Phase 
Shift Across the Array [Ref. 6:p. 43] 
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Figure 10. Steered Directivity Pattern [Ref. 6:p. 44] 
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A planar array has the advantage of resolving the 
azimuthal and the elevation angles of arrival of an incident 
plane wave [Refs. 5,9]. Consider a rectangular-shaped 
planar array as shown in Figure ll; sensor elements are 
arranged in a rectangular grid in the x-y plane. The center 
of the array is usually chosen as the coordinate origin. 

The entire array has M elements in the x-direction with uni- 


form spacing d.. and N elements in the y-direction with uni- 


form spacing d.- The elements are assumed to be point sources. 


The phasor sum of the entire array can be written as: 


jmyp jnyp 
s(t) = }) } y(t) e oc Y (2 Sis 
mn 
where 
dy 
Vy, = 2m (—) sin@ cos¢ (27 
dy 
Us = 27 (—-) Sin8 sind (2.14) 


The directivity pattern of the planar array is given by: 


jmp. jny 
A(f£,0,6)" = Jee e Y = A (£,0,0)A_ (£,8,¢) (25 ESn) 
nen a 


It follows from the model given by equation (2.12) that the 


planar array beam pattern is the product of the array factors 
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of two linear arrays. However, separability of the two- 
dimensional beam pattern is not necessary to ensure the 
proper operation of the LMS adaptive algorithm for a planar 
array. Beamsteering is accomplished by applying appropriate 
linear phase shifts for the row and column elements. For 
the rectangular array case under investigation here, it is 
more convenient to transform the elevation-azimuth (6,6) 
Space to a rectilinear coordinate space (u,v) by the 


transformation: 
u = sin 6 cos $9 (2356) 


V = sin 8 sin o (27) 


The parameters u and v are the direction cosines with 
respect to the x and y axes, respectively. Figure 12 shows 
the alternate diagrams for presenting two-dimensional array 
beam patterns [Ref. 6]. The ranges of the spherical angles 
Poe woe /2 ana 0 < 9 < 27 whereas the ranges of the 


fal einoan Coordinate system are -l1 < u< land-l<v«<is. 


B. THE PLANE WAVE MODEL 
L. Far=Field Condition 
The LMS adaptive filter designed in this thesis 
will provide spatial resolution for a planar array for an 
incident plane-wave field. The plane wave assumption is 


justifiable for a radiating source located in the far-field 


Ii2) 


Pir, 9, 9) 
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Projection of 
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a x-y piane 


td / 
Sensor 
eiements 
Figure ll. Sensor Element Arrangement of a 


Rectangular Planar Array [Ref. 6:p. 45] 





Figure 12. Transformation from Spherical Coorae@nace 
Space to Direction Cosine Space [Ref. 6:p. 46] 
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[Refts. 5,10], due to wavefront expansion. The far-field 


memige tOr a planar array is given by [Ref.,5]: 


2 
m2 (2.18) 





where: 


i ees 72 | 
rn (2) *} (2.19) 


represents the maximum radial extent of the transducer array 
and Ly and hb, are the dimensions of the planar array in 
the x and y directions, respectively. 
2. Propagation of a Plane Wave from a Far-Field Source 
The plane wave solution of the Helmholtz wave equa- 
tion has the form: 


aca) Owe a) 


y (t,x) (ZaZo) 


enemy ats Called the velocity potential, f is the 


PeeqleWeyyie 1S the propagation vector, and r the position 


vector. In rectangular coordinates, 

k = kK % + sy! + kz Zee ) 
and 

Bee xXx yy + 22 CZ APD 


Sk 


For a planar array located at some reference location z = 0, 
the velocity potential is: 


jl2uft+ (k,x+k y)] 
pecc, 2) a Vy (exe = Ae Y (27243) 


For a planar array with discrete sensor elements located 
uniformly in the x-y plane with spacings d.. and d., respec- 
tively, the continuous space variable x can be replaced by 
md, and y replaced by nd,,. If the time signal is digitized 
for computer processing, then the time variable t can be 


replaced by QT. where 2 is the discrete time index and TS 


the sampling interval. To summarize: 
Co ee (2.255 
x > md, (2. 255 
y > ndy (2.25c) 


The corresponding discrete time, sampled space signal is 


given by: 
J27e2T 3 Gnd.,k popes) 
_ S xX Xie oy. 
y (2T_,md,,nd,,) = Ae e (2 lee 
where y(2£T_,md,,nd,) is usually shorted to y(2%,m,n). The 
propagation vector k at z = 0 has two components ise and 


a2 





that can be related to the direction cosines u and v via 


eet 5): 
ae 
kK ys (27a) 
and 
er = 4mV (2527) 
y r 
ieee eueinog Cqudtiom. (2.27) into equation (2.26) yields: 
5 AS , 
j2r£2T. J— (umd, +vnd,) 
y(2,m,n) = Ae e (228 } 
nec 
T Zaeec al, 
y(2\= Ae 2 E29 |} 
represent the time dependence of the signal. Equation (2.28) 
then becomes: 
al 
J (umd, +vnd_) 
y(2,m,n) = y()e (278.0) 


From equation (2.30), it is easy to see that the signal at 
each element location (m,n) has the same time dependence 


but has a different phase due to the element location and 
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the direction cosines associated with the incident angle of 
the plane wave upon the array. The exponential relationship 
of the phase in equation (2.30) also suggests that the 


phases in the x and y directions are separable. 


C. FREQUENCY DOMAIN LMS ADAPTIVE FILTER FOR SPATIAL RESODRI 

The objective of the LMS adaptive filter used in this 
thesis 1s to phase align: the signals from all sensor elements 
such that they add up coherently to realize the full array 
gain. Figure 13 shows a general cophasing scheme for linear 
arrays. Cophasing or phase alignment is done in the frequency 
domain by multiplying the frequency spectrum at each element 
by she proper phase weight in order to cancel out the phase 
due to element location. This is equivalent to a phasor 
rotation in the complex plane. The amount of rotation needed 
to align each sensor element is proportional to the frequency | 
of interest, the direction cosines u, v and the locationwes 
that element. 

1. Phase Weights for the Planar Array 

The total array output 1S maximized if all elements 

in the array are phase aligned. If we let cd(m,n) be the 
proper phase weight at location (m,n), then the phase 


weighted total array output is: 
s(2) = J) ) cd(m,n)y(2,m,n) ga 
ee 3G! 


Substituting equation (2.30) into equation (2.31) gives: 
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Figure 13. Adaptive Cophasing Scheme for a Linear Array 
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OF 
j—(umd,+vnd,,) 
s(2) = y(2)o) ) Cdinpc ‘. a 


mM fn 


(2s 


iit 
2 
= | a (umd, +vnd,, ) 
Cantina ) = e (2.349) 


then the quantity inside the summations becomes unity and: 


Aa = y Cee a) (2 gees 
mn 


Or 


Ss (2) MN y(2) ( 285) 

where s(2) is the sum over all elements and equation (2.35) 
is the maximum signal level possible. .This maximum level is 
achieved by tuning M xN adaptive weights cd(m,n) to conform 


to equation (2.33). The same phase weighting procedure is 





also true in the frequency domain, in fact, the implementation 
of phase weighting is inherently a frequency domain opera- 
tion. Since phase weight equation (2.33) is a function 
wavelength \ and i = c/f, the proper phase weight for co- 
phasing at each element is a function of frequency fe. Thus 


for each valid frequency component in the signal, a different 
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set of phase weights {cd(m,n)} must be generated. Consider 
the DFT with respect to time of equation (2.32): 

2 Gmne +vnd,.) 

IX ‘aan 


S(q) = Y(q) } ) cd(m,n)e (2.36) 
mM fn 


where q is the frequency index and Y(q) is the DFT of y(). 


Tf equation (2.. 33) holds, then: 


ehcp ee cm) (1) = MN Y¥(q) E37) 
Nn n 


Only valid spectral lines will be processed. 
2. The Frequency Domain LMS Adaptive Filter 
The general frequency domain LMS adaptive algorithm 
is derived in Appendix A. Suppose that the time signal z(Q) 
is the reference or desired signal and z. (2) is the normalized 
sum of all signals in the planar array. In the frequency 


th 


domain, the reference signal in the gq DET sbaen) is: 


eZ 1 


Mill “I= *%49 
Ziq) = ) az(R)e (23.8 ) 
=0 
and the estimated output in the frequency domain is: 
a L-l A -j22q 
Z.(q) = } 2a, (%)e 2.38) 
; 2=0 


Bef 


A ie el -j<Teq 
Z.(q) = 2} ge?) } ed, (m,n)y(2,m,n)e (2.40) 
L=0 ig 1 
where: 
§ 1s the time indéex-3))— 0-1 
- q 1s the DFT bin index Go =s0) 12 ol 


1 ais the complex phase weight iteration number. 


The DFT operation with respect to time in equation (2.40) 


can be performed first to yield: 


- eee 
Z. (q) = Sa > cd. (m,n) ¥(q,m,n) ( 2 Sale 
mn 
Wie we. 
b= -jTeq 
Y(q,m,n) = >} y(2,m,n)e (242) 
L=0 


Based on the complex LMS algorithm [Ref. 1], the adaptive 
filter output in the qth bin is given by equation (2.41). 
The error signal is generated by comparing the desired 
(reference) signal to the adaptive filter output (estimate). 


The error is denoted by eu, where: 
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“~~ 


a = Z(q) - @, (q) (27435) 


The estimate, equation (2.41), is formed by applying the 
phase weights {cd, (m,n) } of the ve iteration to each element 
in the planar array. The complex weight cd. (m,n) is updated 


recursively as follows: 


= * 
cd. ,, (m,n) cd. (m,n) + 2u.e.¥ MG Mp tl) (244 ) 
where: 
ig == 90) elie »M-1L 
me = Oeil]. .,N-L 
(*) denotes complex conjugate 
Hy, = feedback coefficient, a parameter that con- 


trols the rate of convergence, algorithm 
noise, and the stability of the algorithm 
[Ref. 4]. ) 

From equation (2.44), it can be seen that the gee weight 
cd. ,, (m,n) may have magnitudes larger than unity. MThis 
growth in magnitude is undesirable for the purpose of 
spatial resolution since spatial resolution depends on the 
relative phase between adjacent elements to resolve the 
direction of wave arrival. Thus, a normalization is neces- 
sary to bring equation (2.44) back to unity. This normali- 


Zoeion 1S: 


ele, 


cd m,n) <« cd (m,n)/|cd,,, (m,n) | (2.455 


a! itl 
These updated and normalized phase weights can now be applied 
to equations (2.41), (2.43), and (2.44) in sequence to 
compute the next set of phase weights. This iterative 
process stops when predetermined criteria are met. At that 
point, the set of phase weights {ed. (m,n) } can be used to 
find the direction cosines of the incident plane wave. 
However, the phase angles of the phase weights {cd (m,n) } 

may have been wrapped around an integer multiple of 27. The 
procedure for phase unwrapping is explained in Sections 4 and 
5 of this chapter. 

The feedback coefficient We in equation (2.44) con- 
trols the rate of convergence and the stability of the LMS 
adaptive filter. Robbins and Monro [Ref. 11] showed that 
the adaptive weights {ed, (m,n) } will converge to the optimum 


result if We is allowed to decrease with the iteration index 


i. The precise conditions are [Ref. lé]: 
WH, > 0 
<i 
Ps a 
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A coefficient He that satisfies the above conditions will 
work as long as the signal and noise inputs are truly 
Stationary but will not be satisfactory for a filter operat- 
ing in a slowly varying environment. Widrow's LMS algorithm 


[Ref. 1] uses a constant value of uw satisfying the inequality 


where a ox is the largest eigenvalue of the correlation 

feeris Of the input. Although this matrix»is» typically not 
known a priori, some bound can be set up by examining equation 
(2.44). Tf stationarity can be assumed, it is possible to 
update HW, every iteration in order to obtain the optimum set 
Of phase weights. In Appendix A a simple method of updating 
the feedback coefficient Ws 1s proposed to improve the per- 


formance of the LMS adaptive filter. 


3. Applying the Frequency Domain LMS Adaptive Filter 
eemosr lLanarsAr nay 


Given a rectangular planar array with M xN elements, 
there are several ways to process the signal from the array 
to achieve spatial resolution in both azimuth and elevation. 


Three different ways are considered in this thesis. 


Al 


a. Orthogonal Linear Arrays 
Two-dimensional spatial resolution is possible 
by just considering one linear array in the x-direction and 


one linear array in the y-direction. A total of M+N-1 ele- 


ments out of MxN are used. This scheme is useful when 
processing time is limited. Figure 14 illustrates the j¢ie@ma- 
of the center linear arrays for this scheme. However, any 


two orthogonal linear arrays in the planar array can be 
used. The recursive equations needed to implement this 
algorithm are divided into two sets; one set for the linear 
array in the x-direction and the other for the y-direction. 


In the x-direction the estimate is: 


“a if M-1 
Zz. (q) = Te c. (m) ¥(q,m,n =H) (2 
ie m=0Q 
where: 
aaa constant y-direction index 
c; (m) = unity magnitude phase weight. 


The error is the difference between the reference and the 


estimate. 


e. = Z(q) -Z% (q) - (2s 


The recursive update for phase weights 1s: 
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Figure 14. A Particular Orthogonal Linear Arrays 
Configuration 
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The phase of the new update is: q 


In the y-direction, the procedure is Similar: 


N-1 
st: _ 

Estimate: ae Gi = ) d.(n)¥(q,m =m, ,n) (2.508 

iL n=0 
|p a ley as = 2Z(2z2) - Z@ 2 son 
ey. Z,) y, ‘v ( ) 

. = * = 

Update: ds 4 fn) d. (n) + Ahy ey .* (g,m morn) (2252) 
Normalization: d,,,(n) + do 47 (n)/|4.,, (9) | (2 Say 


The convergence constants Ly and Me are usually set to be 
equal since the statistics in the orthogonal directions of 
a planar array can be assumed to be the same for the obser- 
vation time of most systems. 
b. Two-dimensional Array 
This scheme uses all M xN elements and therefore 


it can realize the full array gain of equation (2.35). iiae 


% 
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phase weights cd(m,n) are not assumed to be separable. The 
eguations for the LMS adaptive filter in the 1 SSMS ey 


domain are: 


a 1 
Estimate: 2.(q) = gy) ) cd, (m,n)¥(q,m,n) (2,54) 
mn 
BieGOr : a ZHCce) = Z. (q) CARS e)) 
Update: cd. 4 (m,n) ~ cd, (m,n) + 2u,e,Y* (q,m,n) 2925 6!) 


Normalization: cd,,,(m,n) < cd (m,n)/|cd.,, (m,n) | (eT) 


i+l 
c. Separable Two-dimensional Array 
As mentioned in the discussion on planar arrays 
and plane waves, the form of a plane wave suggests that the 
phase of a signal at an element (m,n) is separable. This 


scheme then uses the separability property 
cane) = ce (m)d (n) oS) 
to implement a two-dimensional LMS adaptive filter. All 


MxN elements are used but only M+N phase weights need to be 


updated recursively. The equations are: 


Estimate: Zz. (q) = d, (n) ) c, (m)¥(q,m,n) (259) 
Nn m 
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Jcige nage e. = Zig Z. (a) (2-600 





asees C54) (mM) = GaGa) ayes. d,(n)¥(q,m,n))* (2 6a5) 

d5 7 (7) = d. (n) ot 2uiei c.(m)¥(q,m,n) )* (2.62) 

Normalizations: c,_, (m) Cc. 41 (m)/|c 47) | (2 an 
nl 

di4,(m) + 45.) (n)//4,,, () | (2.64) 


All three of the aforementioned schmes are implemented and 





their results compared. At the start of all three algorithms, 
the initial phase weights are set to the boresight of the 
planar array, i.e., magnitude equal to unity and phase angle ) 
to zero. The normalization of phase weights to unity forces 
the spatial transfer function of the LMS adaptive filter to 
have unit magnitude. The steady state phase response is 
designed to phase aati all sensor elements in an element- 
by-element fashion. More discussion on this topic can be 
found in Appendix A. 


4. Extracting Estimates of the Direction Cosines u and 
v from Phase Weights 


To extract u and v from the orthogonal linear arrays 
and the separable two-dimensional array cases discussed in 
Section C.3.a. and C.3.a.c., only M elements in the x-dimegs 


tion and N elements in the y-direction need to be considered. 
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a. Direction Cosine Estimates for Linear Arrays 
Consider that in cases 3a. and 3c., the phase 
weights c(m) , d(n) have reached a steady state. The 
objective at this point is to relate the phase angles of 
these two sets of phase weights to their respective direction 


cosines. Let Ey. (m) be the phase angle of c(m) and Bap be 


the phase angle of d(n), i1.e., 

E _ (m) = tan (SR (2.65a) 
and 

cyim) = tan 2aniginn 2.65) 


It can be seen from equation (2.33) and using the concept 


of separability that: 


e(m) = - Su may) (2. 66a) 
ee ee 
Spy = 7 (v ndy) (2.66b) 


Solving equation (2.66) for u and v yields: 


A ~AE,, (m) | 

u = ~2amdy ’ m a 0 (2.67a) 
x eal 

V = ~2and ’ nN # 0 (eo OD} 
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Thus, in the x-direction, there are M-l estimates of u; 
while in the y-direction, there are N-l estimates of v. 
To find an estimate of the direction cosines from equation 


(2.67), one needs to take an arithmetic average of possible 


estimates: 


is a: 
u = = > —_— (2.68a) 
M-1 = 2mmd., 
z N-l -A&_(n) 
= 1 y 
veer l a (2.68) 
n=] Y 


Equations (2.68a) and (2.68b) will be referred to as the 
point by point method. Another way of finding = and . makes 
use of linear regression [Ref. 13]. Consider equation (2.66) 
where Ey (m) is a linear function of m with slope equal to 
Tad and 2g is a linear function of n with slope equal 
to -“Tyd,,. Using a linear regression fit of M data points 


vs. the element number m, the slope and intercept of the line 


Ey (m) can be calculated. The same procedure can be used 
zone Fey Let the slope in the x-direction be Sy and the 
Slope in the y-direction be So: Then: 
_. 
al = zy udy ; (2.69a) 
s, = = 72a (2.69b) 
2 r Y 
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Thus, Wu and v can be solved by rearranging equation (2.69) 


fornne form : 


A -AS) 
1B} - 2d, (22/70a} 
A “AS5 
V = 27d, (2270b) 


The estimates 5 and v obtained from the linear regression 
method represent the best linear least-squares fit of the 
observed data. 
b. Direction Cosine Estimates for Planar Arrays 

The two-dimensional array discussed in Section 
3.b of this chapter does not require the phase weights to 
be separable. Consider that the set of phase weights 
{cd(m,n)} has reached a steady state. Let Ssey (Men) be the 


phase angle associated with the cd(m,n). Then: 





_ = mlcaum, nm) | 
Sexy (Men) = tan “{eetéatmen)l (2.71) 
From eguation (2.33), it can be seen that 
| Se) En md., + Send ) (ea 2 
| oye r X Y 


In general, Exy Men) can be a more complicated function of m 


and n. For instance, in the near-field problem, one needs to 
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modify equation (2./2) to contain quadratic phase uremic. 
account for wavefront curvature [Ref. 5]. For the plane 
wave model, equation (2.72) adequately describes the phase 
weights needed to steer the directivity pattern of the planar 
array to the direction corresponding to i and = 

The point by point method is applicable here 
TS) se alsovel u and 7 given a steady state phase angle. However, 
a linear regression fit of Seay, PD) vs. the element index 
mand n appears to be more suited to this problem. Rewriting 


equation (2.72) in the form of the equation for a plane in 


three-dimensional space yields: 





2 ui = (=a ud,,)m af (ro vdy)n (27ee 








Equation (2.73) describes a plane with slope -Tua, in the 


m-direction (x-direction and slope - va, in the n-direction 


(y-direction). Again, let the slopes be: 
eee! on a (2.74a) 
il X X 
o  2 =a (2.74b) 
2 7 Y 


Thus, equation (2.74) is identical to equation (2.70) jam 


the direction cosine estimates u and v are found by equation 


(2.71) sen 
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A “AS, 
Vv = 2rd, OLN FE ley) 


This result is not surprising since the exponential represen- 
tation of the plane wave, equation (2.20), 1s inherently 
Marable. 
5.  Unwrapping the Steering Phase Weights 
a. Linear Array Unwrapping 
The proper phase weights for beamsteering are 


given by equation (2.66): . 








So, Ue =2t iy md.) (2. 66a) 
g(n) = “22 (v nay) (2. 66b) 


Consider a 7-element linear array lying in the x-direction 
with the center element as the reference element. The element 
MeGe< sEnene Guns from m=|]= —3, —-2, -l, 0, 1, 2, 3. If u is 


equal to 0.55 and d.. = \/2, then equation (2.67a) reduces to: 


(Zen) 


sy 
3 
iI 
( 
3 
cS 
= 


=O. 95. ih mi 


oll 


The following table shows the required phase weights Ey (m) 


needed to steer the beam to u = Q.55. 


TABLE l 


PHASE WEIGHTS FOR BEAMSTEERING 


m ~3 ~2 -1 0 1 2 3 
Em) 1.657 eal 255i 0 -.557 — Laie —oom 
€,- (m) == Soil = SIN) A) (I) 0 -.557 Oa Fo) 

J&,, (m) 
The phase factors {e } are a set of complex numbers in 


the complex plane. The angle of any complex number must lie 
within a 27 interval. The interval chosen here is [-1T,mT]. 
This means that any angle Ey (m) that is outside the range 
[-7t,17] will be wrapped around to an angle Ey. (m) that is 


inside the range. This property can be shown as follows: 


el (Bt2nk) _ gJB,j¢mk _ 38 (2.76) 


since 


oJ 4mk = eT 


therefore 


alle has modulo 27. 


a2 





Referring back to Table 1, the observed angles 65, (m) are 
wrapped. If no processing is done to unwrap the observed 
angles, the equations derived in Section C.4 of this chapter 
will not apply for all permissible values of a and = For 
small values of u or v, the needed phase weights do not wrap 
around but the spatial range of interest is severely restricted 
[Ref. 4]. In tracking systems, the above restriction in look 
direction can be Pett Ficd Since a crude target direction is 
usually provided by a search array. The maximum spatial 
window of an M element linear array without the unwrapping 

of the steering phase weights can be calculated. For example, 
given that the element spacing is d,. = /f/2, and the maximum 
permissible magnitude of Ey, (m) is tm in the range [-1,1T], then 


equation (2.66a) becomes: 


ce eT = | rum | (Zari 


Or 


_ 


lum | 
max 


For the 7-element linear array described here, the maximum 


value that the index m can have is 3. Therefore, 
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in. ciemer ale, 
(2. vim 


where M is the number of elements in the array. The corres= 
ponding angular coverage in terms of the polar coordinate 


angle @ is (for © =o) 


@ = sin (ul...) (2.79) 
The total angular coverage iS 20. HOr the ease) Ome aaa 
element linear array, this corresponds to a coverage of about 
40° out of a range of 180°. Figure 15 shows the expected 
phase angles required for beamsteering vs. element number. 
Figure 16 shows the wrapped angles vs. element number. It 
should be noted that either set of these angles (phase 
weichts) will steer the beam to the proper direction. The 
difficulty with the wrapped (observed) angles 1s) that ame 
direction cosines u and v cannot be directly estimated using 
the methods developed in Section C.4 of this chapter. In 
order to unwrap the observed angles in Table 1, consider 


equation (2.76). It can be seen that the observed angles 





differ from the angles generated from equation (2.66) by 
an amount of +2nk where k = 0,1,2,.... In order to unwrap 
the observed phase angles, it is necessary to find out which 


elements’ phase angles have been wrapped around and by what ! 
I 


54 





PHASE ANGLE IN RADIANS 


Figure 15. 


ELEMENT INDEX 


Phase Angles for Beamsteering 


fe) 5, 


PHASE ANGLE IN RADIANS 
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Figure 16. Wrapped-Around Phase Angles 
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integer multiple of 27. Let w(m) be the unwrap factor such 


that: 


Pm Fe. (im) =  € 


J . x (M) + w(m) Cow) 


Table 2 shows (for u = 0.55 and d.. = X/2) how the implemen- 
tation of equation (2.80) will yield the phase weights 


required for spatial resolution. 


TABLE 2 


PHASE UNWRAPPING 


m -3 —2 1 0 1 2 3 
Bm) IGS 1 ddr woon Or = OS ab diy 1 le Si) 
E- (m) -.35m -.9T 551 0 -.551 0.9m 25m 
w (m) +27 +27 0 0 0 ~27 —27 
E(m) =: 1.65 L.1n 551 0 -.55t -l.ln  -1.650 


Observe from Table 2 and equation (2.66) that the phase 
angles for elements m = -l and m = 1 do not wrap around as 
long as d.. < x. Recall also that the center element is 
chosen as the reference element. It is possible then after 
a set of steady phase weights has been computed that an 
estimate of the direction cosine can be computed using the 
phase angles at elements m = -l and m ='l in equation (2.67) 


below. 


a 


u (m) sand (25a 
X 
me) 1 
d.. = A725 
A ~&,, (m) 


The estimate of direction cosine u using only the information 
(phase angles) at m = +l is denoted as Hs where 


Geese (ul) Pe i) | (2eos 


Equation (2.82) will yield a good estimate of u as long as 
dy. Ss x. The next step in this process of phase unwrapping 
is to use the result from equation (2.82) to generate a see 


of projected phase weights LE, (m) } using the following 


relation: 
F o(m) = Se ( 2B» 


Recall from equation (2.76) that a phase angle outside the 
rangle [-1,71] is mapped to an angle within [-1t,7t]. By 
examining the magnitude of the projected angle, it is. possible 
to determine how many multiples of 27m were lost due to the 


modulo 2m property of complex numbers. The sign of Ug and 
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the location of the element m determine the sign of the 
unwrap factor w(m). For linear arrays with an odd number of 


elements, the unwrap procedure can be described as follows: 

ie oe om) | <. kt (2.84) 
and if |&(m)| does not lie between the limits in equation 
(2.84), the value of k is decremented by two and the inequality, 
equation (2.84), is tested again until |&(m)| falls within a 
27 interval, then 


fwim) | ° = (k=l) (228)5) 


For an M-element (odd) array, the initial value of k is 


- . The sign of w(m) is determined by: 
Sign w(m) = - sign (m) sign ee (2.86) 
where: 
M=] M=1] 
m = — (x7) pe ee 2,71 ,0,1,.--, (57) 


i.e., using the center element as the reference. A similar 
procedure can be utilized in an array with an even number of 
elements. The reference used should then be the point be- 


Gwecem the two center elements since there is no need for the 
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reference point to coincide with the location of a sensor 
element. This choice of reference takes full advantage of 
the resulting symmetry. From Table 2, it can be seen that 
only half of the elements need to be examined using equations 
(2.84)-(2.86). The unwrap factors {w(m)} for the other half | 
are simply the negative of the first half. This unwrap 
procedure 1s good for the two array configurations discussed 
in Sections 3a. and SCanos this chapter. The unwrap proce— 
dure for the two-dimensional array 1S similar but the compu- 
tation is a little more involved. The unwrapping procedure 
for the linear array in the y-direction is identical wien 
the substitution d.. = ds ieee C ss E,, (m) > ST eter 
b. Two-dimensional Array Unwrapping 
The proper phase weights to steer a beam to 


AN ”~ 


(u,v) are given by: 
Se nil sia = + (au md, + v nd.) , (25375 


The unwrapping procedure is best illustrated by considering 


nN MA 


the following example. Consider the case (u,v) = (-0.7,+0.7) 
and d. = d,, = \/f/2. Equation (2.87) then becomes: 
Sy = -7T(um + vn) (2.88) 


= = 0. 7 ie ee eine 
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Assume also that we are given a planar array with M elements 
in the x-direction and N elements in the y-direction (the 
corresponding element indices are (x,y) +> (m,n), with both 
fmena N equal to 53. For the sake of symmetry, let the 


element indices run as follows: 


iteetcnis Case, for M= 5, 


pelariy, for the orthogonal direction, 
igh = I 0) Lae 


The obvious choice of reference element is (m,n) = (0,0), 
i.e., the center element. The desired phase weights for 
this example are given in Table 3. 

The phase weights marked within the two triangles 
in Table 3 will be wrapped around, so the actual observed 
phases when the phase weights reach steady state are ‘tabu- 


ieiwead in Table 4. 


LL 


DESIRED PHASE WEIGHTS {&xy(m,n)} GIVEN 


TABLE 3 


THAT u = -0.7,v = +0./7 


-2.8T7 -2.17 -1.47 -O.77 
= 2 loa -1.47 -QO.7T7 0 
-1.47 -QO.7T7 0 Oe 7 
-O.7T7 0 O.227F dle, CS 
0 Oe7 i ie ob 2 «em 
-2 —] 0 i. 
TABLE 4 


OBSERVED PHASES OF THE STEADY STATE 
ADAPTIVE WEIGHTS Exy (Mn) 


= eM) 10 di 0.67 =O) rene ait 

- lkg 0.67 = aenin 0 

5 (eyith =O al 0 Ocal 

rin 0 0) es 7 Oars 
Oren =P Ol Oren 
—l 0 ” 
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Let w(m,n) be the two-dimensional unwrap factor, 


unwrapped phase is: 


eatin, 1) - 9 = Sy WEY = 


E q 


i) FW) 


KY 


the 


(2210.9) 


By comparing CE (msn) J irae 3s with fey (men) 3 in Table 


4, it can be seen that w(m,n) 


lated values in Table 5. 


must be equal to the tabu- 


TABLE 5 


UNWRAP FACTORS 


n 
Z = Ait = 21 
a = i = 27 
0 = 2H 0 
—ih 0 0 
=2 0 0 


To generate w(m,n) , estimates of 


u and v must be computed using elements at m = +1 and 


-27 


25 


w (m,n) 


27 


21 


2 


et 


EHeudil rection cosines 


n= tl. Let ae and We be those estimates obtained using the 


point by point method below: 


6 3 


N _ L “~~ “~~ 

a= xfu(l) + u(-1)] (2.90) 

vy. = S{v(1) + v(-1)] 2.91 

es = 57 V (2.905 
ilove Gl = d., = }/2, equation (2.88) can be applied to com- 


pute the projected phases te. (mon) by the relation: 


z (m,n = -T u m + v n} fon allieimeands mn 2,92 
Eyey (MD) {u, : (23s 
The set of projected phases (angles) are then examined to 
decide the proper unwrap factor for a particule element 

(m,n). The logic is as follows: for each n, all elements m 
are examined. 


“Ww 


Check (k-2)7 < |& (m,n)| < kn (2.93) 


xY 


and if | Ess (myn) | does not lie between the limits, then kus 


decremented by two and equation (2.93) is tested again. igs 


nN 


Ey (myn) | does lie between the limits, then 
Woes, (msn) | = Sas 7 (2.94) 


This procedure applies when M is odd and the initial value 


of k is (—3—). The sign of w(m,n) is determined by: 
sign[w(m,n)] = ~sign[sign (m) sing (ug) +sign (n) sign(v,)] ( 2a) 


64 





where: 


M-1 MIS 
m — erm oe 07 l, =. , (>) 


a= oe iyo 1,..., (—>—) 


Pwieecne Next value Of n, equations (2.93), (2.94), and (2.95) 
are repeated for all values of m. This continues until the 
last value of n is reached. It is possible though to examine 
Pmeenalt Of the planar array since the unwrap factors for 
the other half are the negative of that of the first half. 
To ensure symmetry about :We reference element, it is neces- 
PelmecO rotate the phase angle at the reference element to 
zero. This can be accomplished by multiplying the phase 
weights of all M WN elements by the complex conjugate of 
the reference phase weight, l.e., 

cd. (m,n) cs cd, (m,n) cd, (mj .n,) Soar ace 
where m and n are the indices locating the reference element. 
This operation will ensure that the unwrapping procedure 
will work properly. For linear array phase weights, this 
phase centering should also be completed before unwrapping. 


The equations are: 


* 
era) ee c. (m) c. (m,) (2.97) 


53) 


d. (ni) = GG tee tae (2.98) 
where mo is the index locating the reference element of the 
linear array along the x-direction and oe is the index 
locating the reference element of the linear array along the 
y-direction. 

6. Summary 

The original complex LMS adaptive filter [Ref. 1] 

requires three necessary modifications to make it useful Gen 
estimating the direction cosines of an incident plane wave. 
Without the following modifications, accurate spatial resolu- 
tion is not possible: 


- normalization of the adaptive complex weights (phasors) 
to unity magnitude after each iteration. 


- unwrapping of the observed steady state phase angles 
to extend the spatial coverage to the full range of 
Wi naniGhs siz 

- allowing the feedback coefficient uw to decrease with 
increasing iteration number to achieve convergence to 
an optimal set of phase weights and to realize a 
Teese Eaileen. 

D. = NOLSE MERE: 

In the SONAR environment, the ambient noise field is a 
composite of many different noise sources. Therefore, using 
the Central Limit Theorem [Ref. 13], the ambient noise can 
be modeled quite adequately as additive white Gaussian noise. 


Intentional jamming is not considered in this thesis. How- 


ever, the use of an adaptive filter to place a null at the 
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Spatial location of a jamming signal has been studied exten- 
sively [Refs. 2,4,15]. The noise corrupted received signal 


1s given by: 
mae sip ee my) + n(2,m,n) 299) 


where y(2,m,n) is the sampled signal represented by equation 
(2.26) and n(2,m,n) is white Gaussian noise time samples 
with zero mean and noise power a 


amplitude (see equation (2.26)), then the signal-to-noise 


If A is the signal 


ratio is given by: 


ne 
on 
Or 
(SNR), = 10 log,)(SNR) 4B (2.101) 


The performance of the frequency domain LMS adaptive filter 
was tested for various signal-to-noise ratios. The noise 
environment is assumed to be stationary during the look 


interval of the SONAR system. 
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ITI. LOW FREQUENCY PASSIVE SONAR TARGET TOCA ZAT Te) 


The objectives in passive SONAR are to detect and possibly 
classify noise sources in the ocean. Most of the noise gener= 
ated by vessels is concentrated in the low frequency range 
(30-1K Hz). Acoustic energy in this frequency range is 
capable of long range propagation and, as a result, most 
long range detection systems in the underwater environment 
Operate at these low frequencies. Propeller cavitation, 
machinery noise, and wake are the major sources of such noise. 
Under certain situations, very long range detection has been 
demonstrated in this frequency range. Ina long range 
detection scneario, Knowledge of the elevation/depression 
angle is necessary to resolve potential range ambiguities 
in convergence zone (CZ) problems. Therefore, the use of a 
planar array is well justified. 

A computer program was implemented using the VS APL 
language to test the performance of the LMS adaptive filter 
in a noisy passive environment. The array size used in the 
Simulation is small compared to most modern systems but the 
other parameters are set to simulate a very realistic pas- 
Sive SONAR. The APL language is used because of its large 
library of advanced signal processing functions and its 
interactive mode of operation which allows for rapid program 


development. The major disadvantage of the APL language is 
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= ee 


its slower speed of execution compared to running a compiled 
Fortran program that performs the same functions. This is 
not a serious problem, however, for the research of this 


thesis. 


A. PROBLEM STATEMENT 
A low frequency acoustic signal from a far-field source 
is received by a planar array in a noisy underwater environ- 


Meme. The moise is assumed to be Gaussian and uncorrelated 


with the signal. It is also assumed that the noise between 
elements is uncorrelated. The planar array has a square 
structure with M = 5 elements in the x-direction and N = 5 


Mlenents in the y-direction. The element spacing is set to 
one-half the wavelength of the maximum frequency of the 
system's operating range. The geometry of the problem is 
shown in Figure 17 [Ref. 18]. The system parameters are: 


- Signal: A cos(2nft), where A is the amplitude and 
f is the frequency 


- Integration time: 0.5 seconds (TO) 

= Frequency recolhuUttonca. 2. HZ (1/T, — £,) 

= Niger of samples: 128 = pt (L) 

- Number of sensor elements in the x-direction: 5 (M) 
- Number of sensor elements in the y-direction: 5 (N) 


- Speed of sound: 1500 meters/second (CO) 
- Sampling rate: 256 samples/second (£.) 
- Frequency range: 0-128 Hz 


- Number of frequency bins: 128 


oe 


TRANS OWGi 


/* 
=< = 
*» 


y 





Figure l/. System Ge@niemn, seo weece ore 
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- Elements spacing: 7.5 meters (d, and d., equal) 


xX 
- Noise model: additive white Gaussian ~ N (0,04) 
- Number of noise samples per sample function: 3200, 


iene Mere = bZo 8x 5) x"5 


The signal used in the simulation is a 100 Hz sinusoid 
corrupted by white Gaussian noise. A total of 128 samples 
are taken for each processing period of 0.5 seconds. This 
corresponds to 256 samples per second which satisfies the 
Nyquist sampling theorem [Refs. 12,16,17]. The maximum 
observed frequency in this case is 128 Hz. The 100 Hz 
Signal will center in bin number 50 for this simulation. 
Peegqueney bin numbers 64-127 correspond to the negative 
frequencies [Ref. 16]. The element spacing of 7.5 meters 
is the maximum allowable separation for the 100 Hz signal. 
to avoid grating lobes in the far-field directivity beam 
pattern. The speed of sound is the speed in the proximity 


O£ the planar array. 


B. SIMULATION 

Given the system parameters stated in Section A, and a 
PZoepoine DEE, a 100 HZ sinusoidal signal will be centered in 
frequency bin number 50. The logical flow graph of the simu- 
lation program is shown in Figure 18. The principle of 
Superposition allows different frequency bins to be processed 
independently. However, if the same frequency is emitted 
from two or more spatial locations, the LMS adaptive filter 


will lock on to the one closest to boresight. Complete 
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Figure 18. Logical Flow Graph of the Simulation 
Program for the Low Frequency Passive 
Detection Problem 


i 





documentation of the simulation program for this case can 
pewEounce in Appendm@wss. Thesoutputs of the simulation pro- 
gram are the estimated direction cosines u and v. The RMS 


Error is defined as: 


Wa 


a. = ies ~~ hear (Sl) 
where: 
Au = 5 ao (Sie) 
Av = _ = (3i3) 
where u and v are the actual direction cosines. This measure 


of error is consistent with the least-squares criterion used 


in estimating u and v. If the estimate of the spherical 


coordinates (6,6) are required, the following transformations 


wn OMAN 


eee transform (u,v) to (6,6) Rei. lo les 


5 WAZ 


eS) eigen (ne = (Ae (3.4) 


and 


§ (u,v) tan +(v/u) ie) 


nN UA 


Leoniould be motea that the transformation from (u,v) to 


Ce pelommommineat. & particular value of the RMS error 


ree 


in equation (3.1) can be due to infinitely many different 
values of 6 and i that 1s, given an RMS error ¢«, there is 
no unique set of spherical angle estimates. The spherical 
angle estimates (6,6) depend on the values of both ' and an 
Evaluating the total differential of equation (3.4) and 


equation (3.5) yields: 


~ & “w aN 


sic _ udu + vdv ( 3euen 


\careey Gieeeigiets 5 


and 


ee v_ du (355 


Replacing dé by A6, dod by Ad, du by Au and dv by Av results 


in the following: 


NN ras 


ae 2 uAu + vdAv . (3.8) 
\/ (uz4v~) (1 — (uty 1) 
- u | va 
Ad = xp agi4v - 7) (i 
ue+v u 


Cx RESULTS 
Three versions of the frequency domain LMS adaptive 
filter were discussed in Chapter II. The performance of 


each will be presented here. 
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1. Orthogonal Linear Arrays 
lifecomeoherciearloneuces Only MtN-1 elements. The 


parameters used for the simulation are as follows: 


Signal - cos 2m (100)t (bin number 50) 

Number of iterations (I) - 75 

Initial feedback coefficient (uy) - 0.0005 

Initial feedback coefficient why) - 0.0005 

Scale factor for nu, and Ee) - 0.909 (see Appendix A) 


Input SNR at single array element, i.e., SNR at FFT 
input w- O dB 


Incident angles (8,¢) - (55°,35°) (see Figure 17) 


The corresponding direction cosines (u,v) - (0.67101,0.46985) 


megures 19 and 20 show the convergence characteristics of 
the complex LMS adaptive filter vs. iteration number 1. The 
solid horizontal lines are the true values of direction 
cosines u and v. The oscillations in the beginning of the 
adaptation are due to the large initial values of be and Uy: 
The feedback coefficients a and Uy are scaled down recur- 


sively by the scale factor via the rule described in Appendix 


A which is rewritten below: 
7 a) = wy (i-1)o | (3104 ) 


and 
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Arrays Configuration in Estimating the 
Direction Cosine u at an Input SNR of O9aGB 


Convergence Characteristics of the LMS 


PiGupeerl oe 


[ae 


ee ee ent enon fenton ten ton tension 





=e oa a =~ awe & == _-_ = =» — = aw oe == _—=_ =o ay CE el =o SS oe == = = 


J Ao ee) i ne a 


5 apa 


-_==_ = = =o eo —=E eo ay awe ee oe = == @ a= map = a= a = a a a= == = 2 & == 












i - =a = eel 





pn eet nas et a 





8r'°0 9+°0 vr'O crO0 
A ANISOO NOILOSYIO 








0 


4 
ITERATION NUMBER 


Adaptive Filter using the Orthogonal Linear 
77 


Arrays Configuration in Estimating the 
Duneeceronecosine vy at am Input SNR of 0 dB 


Convergence Characteristics of the LMS 


Figure 20. 


Lot ye epee (3. iy 


Convergence to the steady state occurred in fewer than 35 
iterations in both cases. Figure 21 shows the error in 
estimating u and v using equation (3.1) vs. iteration number. 
The use of a constant feedback coefficient requires a much 
longer iteration period and results in less accuracy (see 
Appendix A). The chosen initial values of Ly and Me, are 
outside the bound of the convergence coefficient described 
in Appendix A. However, the convergence coefficients are 
decreased rapidly using equation (3.10). This choice of 

Hy and iy shows that the LMS adaptive filter with decreasing 


convergence coefficient(s) is very robust. Table 6 summarizes 


the simulation results for a particular sample function of 


noise. 
2. Two-dimensional Array 
This version uses all M xN elements and the complex 
weights are not assumed to be separable. The caraaeee 


used in the simulation are: 


Signal - cos 27(lLOQRe  (biteaember 50) 
Number of iterations (I) - 75 

Initial convergence coefficient - 90.001 
Scale factor (a _) - 0.909 


Input SNR at single array element, i.e., SNR at FFT 
Input = Opas 
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TABLE 6 
SUMMARY OF COMPLEX LMS ADAPTIVE FILTER 


PERFORMANCE FOR THE ORTHOGONAL LINEAR 
ARRAY CONFIGURATION AT SNR = OQ dB 


Simulation Values 





u 0.671010 
V 0.469846 
: | O6S56 75 
+ 0.466038 
Ausu-u 2.4665 x10° 
Av=v-v Bei «10 = 
ce = (Au*t+Av~) - 0.024957 
s) DS a 
: 54.783° 
ob 55 
$ 357 ols 
Incident angle (6.0) 5]. oo) 
Corresponding (U,V) gee Ocoee Or 0 =n e2) 


Figures 22 and 23 show the convergence characteris- 
tics of the two-dimensional array configuration of the 
complex LMS adaptive filter vs. iteration number 1. 

Figure 24 plots the RMS error in estimating u ama 
v (see equation (3.1)). ‘The results of the simulationseem 
a particular sample function of noise are summarized in 


Table 7. 
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TABLE 7 


SUMMARY OF COMPLEX LMS ADAPTIVE FILTER PERFORMANCE 
FOR THE TWO-DIMENSIONAL ARRAY AT SNR = O dB 


Simulation Values 


u 0.67101 

V 0.46985 

u 0.68442 

v 0.46569 
Au beac elo 
Av ASG ie 

; 1.404 x10? 

: 55° 

9 55.876° 

$ 35° 

‘ 34.232° 


3. Two-dimensional Array wiletn Separabl Eeaietomese 


This configuration assumes that the complex weights 
are separable, i.e., cd(m,n) = c(m)d(n). The simulation 
parameters are the same as those used in the two-dimensional 
array case discussed in the previous section. Figures 25 
and 26 show the convergence characteristics of this configura- | 
tion. Figure 27 shows the RMS error versus number of 
iterations. The summary for this run is in Table 8. - An 
alternate way to illustrate the performance of the complex 


LMS adaptive filter is a plot of RMS error vs. input SNR at 
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TABLE 8 
SUMMARY OF COMPLEX LMS ADAPTIVE FILTER PERFORMANCE 


FOR THE TWO-DIMENSIONAL ARRAY WITH SEPARABLE 
WEIGHTS AT SNR = O dB 


Simulation Values 


u O67 501 

V 0.46985 

. 0.66304 

: 0 447355 
Au Eee m0 
Av 3.68 x10F- 

€ | 2678 x ome 

Q SS 

g 54.5652 

0) 355 

$ 35.5348 

a Single array element (SNR at FFT input). Figure 28 shows 


these curves for all three configurations. 


D. SUMMARY 

The algorithm that assumes separable complex weights 
shows better performance over the prescribed SNR range. The 
improvement of performance of all three algorithms with 
increasing SNR is evident from Figure 28. Since each noise 
sample function has a total of 3200 independent samples, the 


RMS error versus SNR curve is rather smooth. If several sample 
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functions are averaged, the slight 'humps' in Figure 28 
can be smoothed out further. The RMS error quantity is 


analogous to that of a ‘miss distance’ on a rectangular grid. 


90 





if eer UeormeOMMUNTICATION [Ret. 18] 


The simulation results of two cases are presented in 
this chapter. They are: 

- Homogeneous medium case, in which the output electri- 

cal signal data at each array element is produced by 
the ocean communication channel simulation computer 

mececuamuherelo). Figure 29 shows the system: geometry 
of this case [Ref. 18]. Note that the ray path from 
transmit to receive array is a straight path. 

- Inhomogeneous medium case, in which the ray path is 
bent due to the variable sound-speed profile. Thus, 
the apparent direction of arrival viewed from the 
receive array is different from the previous case.” 
Figure 30 shows the system geometry of the inhomo- 
-geneous medium case [Ref. 18]. 

meen the analysis of Chapter II, it can be seen that the 
LMS adaptive filter should be able to phase align a planar 
Miatay tO ooant in the direction ofarrival. The signals 
used for this simulation are a CW pulse and a LFM pulse. 
These are very common signals used in the SONAR environment. 
Information is carried by the modulation of these pulses. 
The Fortran program used to simulate this problem is docu- 
mented in Appendix C. Blount [Ref. 18] studied the effect 
of model-based cophasing on the probability of detection of 
a single pulse. The amount of cophasing is determined by 
the system geometry and deterministic ray bending. It was 
shown that by applying the phase weights generated by con- 


sidering those factors, the performance of a correlator 


receiver was improved markedly. Analysis done on those 


Ol. 
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oS 


steering phase weights showed that the main beam of the array 
directivity pattern was indeed steered to the direction of 


actual arrival instead of line of sight. 


A. TRANSMIT WAVEFORMS [Ref. 18] 
Two types of waveforms were used £0 test the LMS aaa. 
tive algorithm. 
l. Rectangular-envelope CW Pulse 
The signal presented to the processor is a quadra- 
ture demodulated complex envelope of the CW pulse [Refs. 5, 


ds ie 


~ K jk (27£)t 
wie) ) Ze (4.1) 


where "." denotes complex envelope. 

The pulse repetition frequency is the same as the 
fundamental frequency Ey of the finite (K harmonics) fre- 
quency spectrum from which the pulse is ysnthesized. The 
pulse duty cycle is arbitrarily set to 0.5. The compiles 
Fourier series coefficients Zz used to synthesize the com- 
plex envelope of the CW pulse are obtained from a closed-form 
expression for the complex-valued continuous spectrum. The 
Fourier coefficients are obtained by evaluating the closed- 
form expression for the continuous spectrum at discrete 


frequencies corresponding to integer multiples of the 


94 


SS — eer Eee =—hl(ilt—“‘—SS 


Mieriieved! ame gueney lie tollowing Specific transmit 
Signal parameters were used in all CW pulse simulations: 
- Amplitude (A): 40.0 
ey Duty Cycle=(D): O.5 


- Fundamental Frequency (£0): 20/0 sz 


- Harmonic Values: Zz, = 20 emp [50°] 
Z_1 = 2, = 12.3324 exp [30°] 
Z_5 = 25 = ei EXD FaOeal 
Z_3 = 2, = 4.244134 exp [j180°] 
Z_4 = 2, = 0.000 exp [a 0°-] 
Z5 = Z5 = 2.546479 exp [70°] 


2. Rectangular-envelope LFM Pulse [Ref. 18] 

The complex Fourier coefficients used to synthesize 
the LFM pulse are found using a procedure similar to that 
used for the CW pulse except the closed form expression for 
the complex-valued continuous spectrum of the LFM pulse was 
found by using the method of stationary phase. Officer 
[Ref. 20] describes the method of stationary phase as does 
Papoulis [Ref. 21] who also provides a complete description 
of the LFM waveform. The following transmit signal param- 
eters were used in all LFM pulse simulations: 

- Amplitude (A): 40.0 
eo DUtyeGvebe. (D).: 058 
- Phase Deviation Constant (B): 2356.2 Peaney eis 


- Fundamental Frequency (f): 10 Hz 
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- Harmonic Values: 2Z. = 14.60593 c2oeia 


ee 14. 6.059)3 esc aleel 
Z_5 = 2 = 14.60593 exp [309°] 
Z_3 = 2, = 14.6059 36exp wie 


B. PROBLEM STATEMENT 

A pulsed signal (CW or LFM) is sent to an intended 
receiver in the pect Ts, The signal at the receive array 
has a planar wavefront. The direction of arrival of the 
incident plane wave is determined by applying the frequency 
domain LMS adaptive filter to phase align (cophase) signals 
at all sensor elements an the planar array. The discreue 


time signal at element (m,n) has the form: 


27 


7 J 7 (umd,+vnd,,) 
y (2T_,md,,ndy) = z(L£T_)e (4.2) 
where 
, K jk2nf t 
= 4. 
Z(2T_) 7 Ze (Aas) 
2: time index 


m: element index im the —ditrecc ren 


n: element index in the y-direction 
T: sampling period 
u: direction cosine with respect to the x-axis 
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v: direction cosine with respect to the y-axis 


d: interelement spacing in the x-direction 
d: interelement spacing in the y-direction 
Z: complex Fourier coefficient 

f: fundamental frequency 


K: total number of harmonics 


[Ref. 5] shows that the number of time samples needed to 


completely describe equation (4.3) is L, where: 
I ge Be ae I (4.4) 


mie system parameters are: 


CW LFM 
Integration time (To): 5 ns 100 mS 
Frequency resolution (£0): 200 Hz Om HZ 
Number of samples (L): ik i 
Number of sensors (M): 5 2 
Number of sensors (N): 5 5 
Sampling rate (£2): 2200 samples 70 samples 
ber second per second 

Number of frequency 

bins (Q): lee 7 
Element spacing dyer d.,: 0.1229 meters 0.1229 "meters 
Carrier frequency J aE Solely o Kinz 
Noise model: Additive white Gaussian noise for 


both cases 


Number of complex noise 
samples/pulse: DIS das: 


> 


The carrier frequency EY 1s assumed to be known a priori. It 
is found for the CW pulse that the best direction cosine 
estimates are obtained by processing the spectral line corres- 
ponding to the carrier frequency. This is reasonable since 
the signal-to-noise ratio at that bin is the highest. For 

the LFM pulse, all the harmonic lines still have roughly the 
same magnitude after propagating through the medium. The 
accuracy in estimating the direction cosines is about sities 

same for any harmonic line. Therefore the spectral line 


corresponding to the carrier is processed. 


C. RESURES 
The homogeneous medium case is considered first for both 
CW and LFM pulses, followed by the inhomogeneous case. . 
l. Homogeneous Case 
The parameters for the system geometry (Figure 29) 
are: | 
Speed of sound (Co): 1475 meters/second 
Depth of transmit array (YQ): 1000 meters 
Depth of receive array e) 2500 meters 


Cross range (X -Xo): 500 meters 


>: 3000 meters 





Line of sealehieeeean comic = ie 


True spherical angle 6: 31.81° 


True spherical angle $¢: -108.4° 
Direction cosine u: -0.1666 
Direction cosine v: -0.5000 
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The three array configurations of the frequency domain LMS 
feo et otaiter were then applied to the CW and LFM pulse 
cases to estimate the direction cosines. The averaged re- 
sults for one pulse corrupted by 100 sample functions of 
noise at an input SNR of 0 dB at a single array element for 
the CW pulse, homogeneous case are presented in Table 9. 
The initial convergence coefficient is 0.5 and the scale 
Hactor Ot is 0.909 for both CW and LFM waveforms. 
TABLE 9 
PERFORMANCE OF COMPLEX LMS ADAPTIVE FILTER FOR 
SPATIAL RESOLUTION, 100 ITERATIONS, INPUT 
- SNR = 0 dB FOR CW PULSE, HOMOGENEOUS CASE 


Algorithm 


2-dimensional 
with separable 


orthogonal 


linear arrays 2-dimensional 


phase weights array 
u 0), 1666 =) 16% =O 1666 
v E5000 BEET -0.5000 
u = aly =. 160 = GRe7.0 
v =) ASG -0.5003 = 5082 
Au -0.00048 -0.0006 -0.0004 
Av 0.0243 ~0.0003 -0.0032 
; 0.0247 -0.000624 0.0032 
Q 31.81° 31.81° 31.81° 
-108.4° -108.4° ~-108.4° 
; 30.37° 31.84° 32.02° 
A -109.8° -108.5° -108.4° 
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The rms error versus Signal-to-noise ratio curves for all 
three array configurations are presented in Figure 3l. 

The simulation results for the LFM pulse in the 
homogeneous medium at an input SNR of O dB at a single array 
element are summarized in Table 10. A total of 100 differ- 
ent noise sample functions were used to corrupt the received 
Signal. The tabulated.values in Table 10 are the ensemble 
average values. Figure 32 illustrates the rms error versus 


SNR plot for the LFM pulse in the homogeneous medium. 


TAB ee, eG 


_ PERFORMANCE OF COMPLEX*®EMS ADAP TRV Eee Deere or 
SPATIAL RESOLUTION, 100 ITERATIONS, INPUT 
SNR = 0 dB FOR LFM PULSE, HOMOGENEOUS CASE 


Algorithm 
orthogonal 2-dimensional 2-dimensional 
linear arrays with separable array 


phase weights 


u -~0.1666 ~0.1666 -~0.1666 
Vv ~0.5000 -~0.5000 -0.5000 
u ~0.1096 -0.1611 -0.1462 
v ~0.2695 --0.4096 -0.4079 

Au 0.0570 0.0055 0.0204 

Av 0.2305 0.0904 0.0921 
E 0.2375 0.0905 0.0943 
9 31 73m 31.81° 31.81° 
6 -108.4° -108.4° -108.4° 
6 16.91° 26.11° 25.68° 
4 =a eee Silt, 5° =109 76 
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The performance of the LFM pulse is worse than that 
of the CW pulse since the number of samples per pulse is 
seven where the CW pulse has eleven time samples per pulse. 
Consider also the Fourier coefficients given in Section A of 
this chapter. The magnitude of the spectral line correspond- 
ijegeto the Carrier £or the CW pulse is 1.3/7 times larger than 
that of the corresponding spectral line in the LFM pulse. 
However, for both waveforms, it can be seen from Figures 
31 and 32 that the two-dimensional array with separable 
weights has the best performance, followed by the two-dimensional 
array and orthogonal linear arrays, respectively. For the 
CW pulse, accurate spatial foe can be Obtained ror 
SNR greater than -3 dB. For the LFM pulse, the SNR required 
itemeanoout 3 dB. 

2. Inhomogeneous Case [Ref. 18] 
The parameters for the system geometry (Figure 32) 
aces 
- Speed of sound (Co): 1475 meters per second 
- Gradient: 0.017 per second 
- Depth of transmit array (Yee 1000 meters 
- Depth of receive array CY): 2500 meters 
- Cross range OS -X0): 500 meters 
- Line of sight range |r oe 3000 meters 


- True spherical angle @: 30.10° 


- True spherical angle ¢: -~-109.4° 
- Direction cosine u: -0.1666 
- Direction cosine v: -0.4731 
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Figure 30 shows that the path traveled by the acoustic rays 
1s bent. The LMS adaptive filter is able to Beacon the 
actual direction of arrival but not the true line orecicgne 
direction to the transmit array. The initial convergence 
coefficient u is set equal to 0.5 and the scale factor Oe 

is set equal to 0.909. Both the CW and LFM pulses are con- 
Sidered for each of the three array configurations. Table ll 
Summarizes the performance of the LMS adaptive filter applied 
to the CW pulse case at an input SNR of 0 dB at a single 
array element. Figure 33 shows the decline of rms error as 
the signal-to-noise ratio is increased. All three array 
configurations are included in the plot for comparison. 

The performance for the LFM pulse case is tabulated 
in Table 12 for an input SNR of 0 dB at a single array ele- 
ment. Figure 34 illustrates the performance of all three 
array configurations versus signal-to-noise ratio for the 


LFM pulse case. 


D. SUMMARY 

For both the homogeneous and inhomogeneous medium cases, 
the complex LMS adaptive filter performed as expected. The 
array configuration that assumed separability of the complex 
weights consistently demonstrated better performance than 
that of the orthogonal linear arrays and the two-dimensional 
array. The superior performance can be attributed € the 
fact that equation (2.24) which describes the reception of 


a plane wave by a planar array is separable. Therefore, by 
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assuming separable complex weights, the system is set up to 
match the physics of the problem. The performance of this 
pulse communication system can be enhanced by increasing the 
pulse width, taking more time samples per unit time, uSing 
high resolution spectrum analysis, and enlarging the size 


of the array by adding more sensor elements. 


Ie 


V. CONCLUSIONS AND RECOMMENDATIONS 


= —___ 


The frequency domain LMS adaptive algorithm has been 
shown to perform the function of spatial resolution. The 
number of iterations required (approximately 35) to reaeia 
steady state is found to be much fewer than that of a com- 
parable time domain adaptive filter [Refs. 2,4]. The theees 
modifications made to the original complex LMS adaptive 
filter [Refs. 1,4] are: 


- normalization ©£ the complex wergh Game o mamas, 
magnitude after each update 


~- reduction of the magnitude of the convergence 
coefficient for each iteration 


- unwrapping of the phase weights 
It has been shown that the above three modifications enable 
the frequency domain LMS adaptive filter to be applied to 
a multiple element array, to have a fast convergence rate 
and robustness, and to cover the entire angular range of 
8 and ¢ values. 

In Chapter III, a passive low-frequency signal was 
generated to test the performance of the frequency domain 
LMS adaptive filter in the presence of white Gaussian noise. 
The low frequency problem is significant in the area of long 
range detection and possibly long range control and communi- 
cation for ballistic missile submarines. Some form of 


modulation, whether amplitude or angle, will have to be used 


NPIL, 


Pooaminy ys ene dntOrmation. The slow data rate 1s not a major 
Grawback since only certain predefined emergency operation 
codes will be transmitted. The simulation results indicate 
that for an input SNR at a single array element of greater 
than -3 dB, accurate bearing and depressive angle estimates 
can be obtained. 

In Chapter IV, a pulse communication problem was con- 
Sidered. Two types of pulses, ET Ee CW and 
rectangular-envelope LFM pulses, were used. The possible 
applications at this frequency range (5 KHz) are high reso- 
lution SONAR imaging, active SONAR detection, and communica- 
tion-between submerged vessels. The performance of the 
frequency domain LMS adaptive filter was again tested in the 
presence of white Gaussian noise. The simulation results 
indicate that accurate bearing and depression angle esti- 
mates can be obtained for input SNRs of greater than O dB. 
The result here is slightly worse than that of the passive 
low-frequency case due to the shorter data length used for 
the pulse communication waveforms. For the simulations 
in Chapter IV, the number of time samples collected is the 
minimum reguired to satisfy the Nyquist sampling theorem. 

The following three configurations of the frequency 
domain LMS adaptive filter were considered: 

= orthogonal linear arrays 
- two-dimensional array with separable phase weights 


- two-dimensional array 


ia fal 


In all simulations, the configuration that used the separa- 
bility assumption produced the best results. This is quite 
reasonable since this configuration uses all the output data 
from the available elements and more importantly, its assump- 
tion of separability matches the physics of the plane wave 
Signal since the phases of a plane wave in the orthogonal x 
and y directions are separable. 

In the course of this investigation, some interesting 
topics for further research revealed themselves: 

- exact target localization (range, depth, bearing, 
and depression angle) if accurate environmental 


data can be obtained 


- -more efficient update algorithm for the convergence 
coefficient 


- implementation of more efficient software 
- application of the steady state phase weights generated 
by the frequency domain LMS adaptive filter to Blount's 


[Ref. 18] correlator receiver 


- addition of a noise reduction system before spectral 
estimation 


- modifications to make the system jam-resistant 
- applications of other high resolution spectral analysis 
techniques to produce frequency spectra (such as 


maximum entropy [Ref. 22], autoregressive, maximum 
likelihood, etc.). 


lie 





APPENDIX A 


THE FREQUENCY DOMAIN LMS ADAPTIVE ALGORITHM 


A. DERIVATION 

Consider the phase aligning system shown in Figure 13. 
The performance objective is for the adaptive filter to con- 
verge to a set of phase weights such that the array output 
Signal will match a reference signal. In the frequency 
Gomain, the output of a linear array of M elements at iter- 


ation i can be written as: 


SOLS, (A.1) 


Tf Z(q) is the reference signal, then the error between the 


reference and the array output is: 
e. = Z(q) - 4, (q) | (72) 


where e. is a complex quantity. 


The mean square error is: 


= [2(q) -2,(q)](2(a) -2,(q)]* (A.3) 


ey 
li 


IZ(q) |? -2, (q)2*(q) -2%(q)2(q) + 12; (a) |? (A.4) 


iba 


Substituting equation (A. 1) aGnto equations a ee. Z. (q) 


yields: 


The complex weight vector that will minimize equation (A.5) 
is computed from the following recursive algorithm [Refs. l, 


6,14]: 


* 
Ci41 ~ Si * *H 96g ee) 
where Hs is the convergence coefficient used in the gradient 
method [Ref. 6]. 
For the two-dimensional case discussed in Chapter II, 


SectiGnws-b, mesiace Coad with cd.) to yield: 


ae . 
CO = Sy a 
where: 
e. = Z(q) - & J J) cd. (m,n) ¥(q,m,n) (A.8) 
i WD - my 2 2 cei tm mM, 


If the complex weights {ed. (m,n) } are separable, then 


cd. (m,n) = Se. (i yreeecn 
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and: 


C544 (m) c,(m) + 2ujeg d.(n)Y¥(q,m,n))* (A.9) 
diyz, fm) = dj(n) + 2 eit c. (m)¥(q,m,n))* (A.10) 
mmere : 
e. = Z(q) - ~~ J)c.(m) J 4. (n)¥(q,m,n) G1} 
at ome ne - ; 


Reference 4 studied the application of the frequency domain 
LMS algorithm to a split array. The error signal is 

the difference between the outputs of the two arrays. It 

is possible to extend this idea to an (M xN) planar array 

by generating an error signal e, (m,n) at each element location 


where: 
e; (m,n) = Z(q) - cd. (m,n)¥(q,m,n) (A 2) 


for all m,n. Although this scheme is not studied in this 
thesis, the implementation is rather straightforward. One 
needs only to replace the error e, with e, (m,n) vig ean = Wl ila ea @ y=: 
equations and to rewrite the code to generate the error. 
This scheme, however, does not have the advantage of 


averaging noisy data points from all of the elements in the 


aes 


array to come up with an estimate. Therefore, it can be 
assumed that the noise performance will probably be inferior. 

The magnitude of the reference signal Z(q) is generated 
by averaging all (M xN) magnitude spectra whereas the phase 
of Z(q) is taken to be the phase of the reference element. 
A thresholding is done to determine which frequency bins 
should be processed. 
B. MODIFICATION TO THE LMS ADAPTIVE ALGORITHM TO MAKE IT 

USEFUL FOR SPATIAL RESOLUTION 

Three modifications are incorporated in the basic LMS 
adaptive algorithm. 

tf. Normalization of Recursive Complex Weights 

The components of the recursive complex weights (Ci, 

d. and cd.) must be normalized to have magnitudes equal to 
one. This restriction forces the algorithm to achieve 
minimum error by adjusting the phase weights only. The esti- 
mation of direction cosines depends on the linear relationship 
of the phases of the signal across the array. Figure 35 
shows that, without the normalization performed after each 
complex weight update, the vector sum composing the estimate 
can be made close to that of the reference. The correct 
values of the direction cosines can only be achieved if each 
component in equation (A.1) has the same phase as the refer- 
ence phasor. All the vector sums shown in Figure 35 have 
the same resultant phase and magnitude as the reference but 


the phases of the components are not equal. It can also be 


as: 
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REAL AXIS 


Different Ways of Adding Up Vectors 


to Obtain the Same Resultant 


Figure 35. 


Ay, 


seen that there are infinitely many ways to achieve a small 
error without the normalization. The LMS adaptive filter 
implemented for this thesis has a unity magnitude across 
the array. It is, in fact, a phase-only filter. With tire 
normalization, each component is forced to adapt to the 
phase angle of the reference. This is the only way that 
Spatial resolution (based on the amount of phasor rotation 
of adaptive weights) can be achieved. 

2. The Convergence Coefficient (u) 

This quantity is sometimes referred to as the mee. 
back coefficient, signifying the analogy to control systems 
[Ref. 6]. For nonstationary environments, a constant con- 
vergence coefficient is usually chosen [Refs. 1,4,23]. 
However, if stationarity can be assumed, the performance of 
the LMS adaptive filter can be improved by using a mono- 
tonically decreasing scheme for implementing the convergence 
coefficient [Ref. ll]. If a constant convergence coeffi- 
cient is used, the bound on the choice of value is given 


by [Refs. 12,24]: 


-l 
0 <- ale eee (Aye 
where 
al Jk 
Anax « JTEIR) = ee y R(O) = (1+1)R(0) 
—~ 10 
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and R is the covariance matrix of the received noise cor- 
rupted signal. If u is larger than the maximum determined 
Sreequation (A.13), oscillation will occur [Ref. 25] and 
the LMS adaptive filter will not converge. If u is too 
small, then the rate of convergence is slow. Figures 36, 37 
and 38 show the convergence characteristics for different 
values of u. 

The recursive LMS algorithm is based on the stochas- 
tic approximation technique [Ref. ll]. The choice of yu is 


optimal if the following conditions are met [Ref. 12]: 


° . > 
Us 0 
limu. = QO 
41700 1 
[ee] 
Dy Baa 
i=l 7? 


peporeeneularm Choice of He that meets the above four condi- 
tions is ia i/7i VRete Lee "9 tne effect of the adaptation 
decreases with the number of iterations and ceases completely 
for large i. Simulations show that the above choice of Us 


is better than using a constant u, however, the choice 
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CONVERGENCE COEFFICIENT TOO SMALL 
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Pe o= 2 1 Agee ) 


where 0 < — 1 and Le is the initial value, appears to 
speed the rate of convergence and also achieve a high 
degree of accuracy. Equation (A.14) does not meet the 


oO 


moneda tion L u, = °, but simulations show that the response 
of the Res daptive filter is fast and accurate for 

Oey < a. < 0.95. Further investigations of the convergence 
characteristics are warranted. The current implementation 
of the convergence coefficient results from experimentation 
with various values of Ot to achieve the best estimates of 
direction cosines. With this scheme, the initial value of 


WH, can be set guite liberally since the value of Ws decreases 


geometrically. 


C. PHASE WRAP-AROUND PROBLEM 

The resolution of the phase wrap-around problem in the 
adaptive phase weights was discussed in Chapter [II exten- 
Sively. The scheme to resolve the phase ambiguity depends 
on the proper functioning of the elements adjacent to the 
reference element. In the event that some of those adjacent 
elements are not operational, the unwrapping scheme will 
have degraded performance. However, since the reference 
element can be any element in the array, it is possible EO 
shift the location of the reference element to a region of 
mivemairray that has sufficient adjacent elements that are 


functioning properly. 
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APPENDIX B 
DESCRIPTION OF SIMULATION PROGRAM 
FOR THE PASSIVE DETECTION CASE 
The VS APL application package worksapce (CMS file 
ADAPTIVE VSAPLWS) contains all the functions necessary to 
implement the simulation discussed in Chapter III. The 
general processing flow is as follows: 
- generate time samples of a plane wave signal of 
frequency f incident upon a (M xN) planar array with 
angles of incidence (6,0) 


- add white Gaussian noise for a desired SNR 


- compute the discrete Fourier transform (DFT) of each 
of the MxN time sequences 


- determine the spectral line with the largest magnitude 
and its corresponding frequency bin number 


- apply one of the three frequency domain LMS adaptive 
filters (orthogonal linear arrays, two-dimensional 
array with separable complex weights, or two-dimensional 
array) 

- compute estimates of direction cosines. 


Usage of the functions are described below. 


A. SIG2D 
SIG2D generates planar array signal at each element 


location” (equaeion 272395 = 
Syntax: YN..«. AGI SiG2D AG2 


YN is a L XM XN matrix where L 1s Ehe number Of eee 


samples, M is the number of elements in the x-direction, and 
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N is the number of elements in the y-direction. AGl is the 
spherical angle @ and AG2 is the spherical angle 6. The 
speed of sound (c), the signal amplitude (A), the frequency 
of the signal (F), the interelement spacing (DX and DY), 

the sampling interval (T), the number of elements (M,N) 

and the number of time samples (L) can be changed by editing 


mime function. 


B. NORRAND 
NORRAND generates independent white Gaussian noise 


samples. 
Syntax: NOISE <« K NORRAND N Nl 


K is the number of noise samples desired, N is the mean 
Nl is the variance. The noise array NOISE must be reshaped 
to conform to the shape of the signal generated by SIG2G 


gz : 
NOISE + LM N QONOISE 


The standard deviation On = 1 is necessary toscale a samplefunc- 


tion of noise with zero mean and variance 1 to a desired 
Signal-to-noise ratio. The signal power at each element is 


at /2 and the noise power at each element is o The input 


=e 
signal-to-noise ratio at a single array element is then given 


by 


ile) 





2 2 
A 
SNR = £- = a (Bee 
ZG 
TN N 
Solving for On yields: 
2 1A 
°y = ‘otsiRy? ac 


where A is the amplitude of the signal and SNR is the numeri- 
cal signal-to-noise ratio. Therefore, a noisy signal with 


the desired SNR can be generated as: 


RN < YN + oe x NOISE (B.3) 


where NOISE has zero mean and variance l. 


CC. DIBA Vices 
DFTWRT computes the discrete Fourier transform with 


respect to the time index for each element (equation (2.42)). 
syntax; YK =< @pETWRry 7) 


YK has dimensions L xM xN; the first index L is now the 
total number of frequency bins. RN is the noisy signal 
generated by adding noise to the output of SIG2D. A total 
of (M xN) L-point DFTs are computed using a radix-2 FFT 


algoritiane 


AG 


D. ADPLMS 
ADPLMS computes estimates of the direction cosines using 
the frequency domain LMS algorithm for the orthogonal array 


configuration. 
Syntax: QTGT ADPLMS YK 


OTGT is the frequency bin number where a valid signal has 
been identified and YK is the output of the Pr DFTWRT. 
A reference signal at frequency bin QTGT is generated by 
calling the function REFGEN. Direction cosine estimates 

u and v are computed in two different loops since in general 
the number of elements M is not necessarily equal to N. 
Estimates of both direction cosines are generated for each 
iteration under the names UHAT and VHAT. The phase unwrapping 
is done by calling the function DC1DX for the x-direction 
and DCLDY for the y-direction. The number of iterations 
(ITER) and the initial convergence coefficients (MUX, MUY) 
can be changed by editing the function ADPLMS. The scale 
factor (SCMU) for decreasing the convergence coefficients 


can be changed in the workspace. 


fe ~OrFEMS 
QFLMS computes estimates of the direction cosines 
using the frequency domain LMS algorithm for the two-dimen- 


Sional planar array with separable complex weights. 


Sines ORE OF EMS 6y K 
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QTGT is the frequency bin number and YK is the output of the 


function DFTWRT. A reference signal 1s generated by the 
function REFGEN. Direction cosine estimates u and v are 
computed every iteration of the LMS adaptive loop. They 
are stored in the vectors named UHATOF and VHATOF. The 
phase unwrapping is done by using DC1DX and DC1DY in the x 
and y-directions, respectively. Recall that only M+N com- 
plex weights are updated for this configuration because of 
the separability assumption. The initial convergence 
coefficient (MUQ) and the number of iterations (ITER) can 
be changed by editing the function QFLMS. The scale factor 
is named SCMU and is stored in the variable list in the 


workspace. 


F. ADPLMS2D 
ADPLMS2D computes estimates of the direction cosines 
using the frequency domain LMS adaptive algorithm for a 


two-dimensional planar array. 
Syntax: QTGT ADPLMS2D YK 


QTGT and YK are the same quantities described in the last 
two functions. A reference signal at frequency bin QTGT 

is generated. Direction cosine estimates u and v are com- 
puted for each iteration and stored in vectors named UHAT2D 


and VHAT2D. The phase unwrapping is accomplished by using 
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MemuiMerron UC2ZD. Lhe number Of iterations (ITER) and the 
initial convergence coefficient (MUG) can be changed by 
editing the function ADPLMS2D. The scale factor can be 
Changed by assigning a different value to SCMU in the work- 
Space. The use of a different unwrapping function is 
required since the complex weights are not assumed to be 


separable. 


G. REFGEN 
REFGEN generates a reference signal at a particular 


frequency bin QTGT. 


oy meas. CORRE = OG! REPGEN YK 


COQREF is the reference signal used in the frequency domain 
LMS adaptive filter. The magnitude of CQREF is obtained by 
averaging the magnitudes of all (M xN) frequency spectra 

in the frequency bin QTGT. The phase of COQREF is taken 

to be the phase of the reference element in the QTGT fre- 


gqueney bin. 


Hee DCL Dx 
DC1DX unwraps the phase weights for a linear array in 


the x-direction. 


Svea: UWHAT[i)) © 9 N DEiDY DV 


LWA, 


UHAT[i] is the direction cosine estimate u of the juan 


iteration, M is the number of elements in the x-direction 
and CV is a complex vector of M phase weights that will co- 


phase the incident signal. 


ie DGS) 
DC1DY unwraps the phase weights for a linear array in 


the y-direction. 


Syntax? VHAT (1) "Nebel yee 
VHAT[i] is the direction cosine estimate v of the oo ltée 
ation, N is the number of elements in the y-direction and 
CV is a complex vector of N phase weights that will cophase 


the incident signal. 


J. Seze 
DC2D unwraps the phase weights for a two-dimensional 


alinigely 
Symeax: DCZEmMED 


CD is the two-dimensional phase matrix that will cophase 
the incident plane wave signal. This function is used by 
the function ADPLMS2D. 

An example of using this application package is shown 


as follows: 


138 


Voie ~ i> 5SEG2ZD 35 
eeeiotos 9 3200 NORRAND O 1 
Severson  — P7605 SvoeesNOTSE 
ae okN 9o< YN +eSCALE™ NOISE 


eee S DE TWRI RN 


= OUGT ADPEMS’ YK 


- QTGT ADPLMS2D YK 


pee oiG. OFUMS YK 


Signal generation 
noise generation 
noise generation 
noisy signal 


CLAS fOmi CO i mequenc V7 
domain 


estimate direction 
cosines 


estimate direction 
cosines 


estimate direction 
cosines 


This sequence of statements generates a plane wave signal 


incident upon a 5 x5 planar array with angles of incidence 


6 = 55° and ¢ = 35°. The number of time samples for each 


element is 128. A noise matrix is then generated and added 


to the plane wave signal and the discrete Fourier transform 


with respect to time is taken for each element in the array. 


The three array configurations for the frequency domain LMS 


algorithm are then used sequentially to estimate the 


direction cosines of the incident plane wave. 


donk 


APPENDIX C 
DESCRIPTION OF SIMULATION PROGRAM 
FOR THE PULSE COMMUNICATION CASE 
The VS FORTRAN prgrams were written to implement the 
pulse communication problem discussed in Chapter IV. Two 
separate programs were written utilizing essentially the 
same subprograms. These programs handle the quadrature 
demodulated complex envelope signals generated by Vos'[Ref. 
19] program. The programs are available on user account 


O218P at the Naval Postgraduate School, Monterey, California. 


+ 


A. PROGRAM ADBFM 

This program is compiled using FORTVS and is designed. 
to run under DISSPLA. It requires a storage capacity of 
1 Mbyte. The following sequence of commands should be used 
to run the program. 


DEFINE STORAGE 1 M 


Tt Cis 
- GLOBAL TXTLIB VALTLIB VFORTLIB CMSLIB IMSLSP NONIMSL 
- GLOBAL LOADLIB VFLODLIB 


= FLLEDEF 04 DISK fname DATA 
(fname is the filename of the date file) 


- LOAD ADBFM 
- DISSPLA ADBFM 
When execution begins, the user will be prompted to enter 


the desired values for the necessary parameters. These 


dbo 2 





parameters are noise status, input signal-to-noise ratio in 
dB at a single array element, number of iterations, spec- 
tral line to be processed, convergence coefficient, scale 
factor for the convergence coefficient, and the choice of 
one of the three array configurations. 

PoOmeoomecamnrayscontiguration, plots of the estimates of 
the direction cosines are generated. The plots for magni- 
tude and phase of the difference between the reference 


Signal and the estimate are also generated. 


B. PROGRAM ERVSDB 
This program computes the rms errors for various input 

signal-to-noise ratios. If a plot of rms error versus SNR 
(dB) is not required, this program does not have to be run 
under DISSPLA. The following sequence should be used: 

- DEFINE STORAGE 1 M (if plot is required) 

moot CMS 

SeebebAL TXTLIB VYVALTLIB VFORTLIB CMSLIB IMSLSP NONIMSL 


= FILEDEF 04 DISK fname DATA . 
(fname is the filename of the data file) 


- LOAD ERVSDB 


- DISSPLA ERVSDB (if plot is needed) 
or START * (if plot is not needed) 


When execution begins, the program will prompt the user to 
enter an initial input signal-to-noise ratio in dB at a single 
array element. It will then ask for a dB step size such 

that the next SNR is determined by the current SNR in dB 


plus the dB step size. A total of nine SNRs are allowed. 


ino 


. For example; if the inatial SNR is -l2 dB ame tie cceo oe 
1s 3 dB, then the program will compute rms errors for the 
set of dB levels {-12, -9, -6, -3, 0, 3, 6, 9, 12}. The 
other parameters such as iteration number, initial conver- 
gence coefficient, scale factor for convergence coefficient, 
and the spectral line to be processed are entered when 
prompted. The program will then ask for how many sample 
functions of signal and noise are to be averaged. Simula- 
tion results show that the average of 50 to 100 sample func- 
tions are sufficient to reduce the variance of the direction 
cosine estimates. One of the three array configurations is 
then -chosen by the user to estimate the direction cosines. 
The screen output of this program is ordered pairs of 
rms errors corresponding to a particular input SNE Giese 
for all nine specified SNRs. A plot of rms error versus 
SNR in dB can be generated if desired (provided that the 


program is run under DISSPLA). 
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